SLAC-PUB-8676 
October, 2000 



Nonlinear Sf Method for Beam-Beam Simulation Q 



Yunhai Cai, Alexander W. Chao, Stephan I. Tzenov 
Stanford Linear Accelerator Center, Stanford University, Stanford, CA 94309 

and 

Toshi Tajima 
University of Texas at Austin, Austin, TX 78712 

and 

Lawrence Livermore National Laboratory, Livermore, CA 94551 



Abstract 

We have developed an efficacious algorithm for simulation of the beam-beam in- 
teraction in synchrotron colliders based on the nonlinear Sf method, where Sf 
is the much smaller deviation of the beam distribution from the slowly evolving 
main distribution /o- in the presence of damping and quantum fluctuations of 
synchrotron radiation it has been shown that the slowly evolving part of the 
distribution function satisfies a Fokker-Planck equation. Its solution has been 
obtained in terms of a beam envelope function and an amplitude of the dis- 
tribution, which satisfy a coupled system of ordinary differential equations. A 
numerical algorithm suited for direct code implementation of the evolving dis- 
tributions for both Sf and fo has been developed. Explicit expressions for the 
dynamical weights of macro-particles for Sf as well as an expression for the slowly 
changing f have been obtained. 
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1 Introduction 



The effects of the beam-beam interaction on particle dynamics in a synchrotron collider are 
the key element that determines the performance of the collider such as luminosity - ■ 
In order to accurately understand these effects, it is necessary to incorporate not only the 
overall collisional effects of the beam-beam interaction, but also the collective interaction 
among individual parts of the beam in each beam and its feedback on the beam distribution. 
The particle-in-cell (PIC) approach J|], P has been adopted to address such a study need 

B, i, B- 

Particle-in-cell codes typically use macro-particles to represent the entire distribution 
of particles. In the beam-beam interaction for the PEP-II || (for example), the beams 
consist of 10 10 particles each. Simulating this many particles with the PIC technique is 
computationally prohibitive. With the conventional PIC code 10 10 particles are represented 
by only 10 3 — 10 4 macro-particles allowing simulation of the beam-beam interaction in a 
reasonable computation time. However, the statistical fluctuation level of various quantities 
such as the beam density p in the code is much higher than that of the real beam. The 
fluctuation level Sp goes as approximately 

p Jy 

where N is the number of particles. Therefore, the fluctuation level of the PIC code is about 
10 3 times higher than that of the real beam. Although this probability is not significant for 
beam blowup near resonances, the higher fluctuation level has a large effect on more subtle 
phenomenon such as particle diffusion. The purpose of the 5f algorithm is to facilitate the 
study of subtle effects and has been introduced in flQfl , [ p"T|j , fll2"| . 



The 5f method follows only the fluctuating part of the distribution instead of the entire 
distribution. This is essentially modeling the numerator of the right-hand side of equation 



1.1). So the 10 3 — 10 4 macro particles are used to represent a/10 10 or 10 5 real fluctuation 



particles in PEP-II beams. This is only one or two orders of magnitude beyond the number 
of macro particles. Such a modest gap between the number of macro particles and the real 
fluctuating particles maybe ameliorated by the standard techniques of the PIC approach, 
such as the method of finite-sized macro-particles [[|, ||. 

PIC strong-strong codes use a finite number of particles to represent the Klimontovich 
equation for the microscopic phase space density (MPSD) fll3"| . In the particular case of 
one-dimensional beam-beam interaction, 

%+P^~ (K(s)x - F{x- s))^j- = 0, (1.2) 
where K(s)x is the usual magnetic guiding force and F(x; s) is the beam-beam force 

F( V ) = ^ P ( S ). (1.3) 
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The electric field E x (x) is calculated from the distribution of the particles of the on-coming 
beam and 8 p (s) is the periodic ^-function with a periodicity of the accelerator circumference. 
The distribution function f(x,p; s) is represented by a finite number of macro-particles by 



1 N 

f(x,p; s) = — ^2 S(x - x n (s))5(p - p n (s)), 



(1.4) 



n=l 



where N is the number of macro-particles. 

The strategy of the 5f method is that only the perturbative part of the distribution is 
followed. The total distribution function f(x,p; s) is decomposed into 



where fo(x,p; s) is the steady or slowly varying part of the distribution and 8f(x,p; s) is the 
perturbative part. The key to this method is finding a distribution fo(x,p; s) which is close 
to the total distribution f(x,p;s). The perturbative part 5f(x,p;s) is then small, causes 
only small changes to the distribution, and thus represents only the fluctuation levels. If 
a distribution f (x,p;s) close to the total distribution is not found or found poorly, then 
Sf(x,p; s) represents more than the fluctuation part of the total distribution; defeating the 
purpose of the method. The ideal situation is having an analytic solution for fo(x,p; s). In 
this case any numerical truncation errors which result from the necessary derivatives of this 
function are eliminated. If an analytic solution cannot be found, then a numerical solution 
needs to be found which is close to the total distribution f(x,p;s) and is slowly varying. 
A frequent numerical update of fo(x,p; s) would also defeat the purpose of the 5f method, 
since the PIC technique essentially does this also. 

The beam-beam interaction can lead to beam instabilities that disrupt or severely distort 
the beam or gradual beam spreading. The higher the beam current, and thus the beam- 
beam interaction, the stronger these effects become. Therefore, when one wants to maximize 
the luminosity of a collider, one needs to confront the beam-beam interaction effects. The 
operation of PEP-II, for example, is critically dependent on the beam-beam interaction and 
optimal parameters to minimize the related beam instabilities are under intense study. 

The paper is organized as follows. In the next Section we present a brief formulation of 
the problem of beam-beam interaction in synchrotron colliders. In Section 3 we develop the 
nonlinear Sf method for solving the equation for the microscopic phase space density in the 
presence of random external forces. The equation for the fluctuating part Sf is being derived 
and its solution is found explicitly in terms of dynamical weight functions, prescribed to each 
macro-particle. In Section 4 we solve the Fokker-Planck equation for the averaged slowly 
evolving part of the distribution. We show that the solution is an exponential of a bilinear 
form in coordinates and momenta with coefficients that can be regarded as generalized 
Cour ant- Snyder parameters. In Section 5 we outline numerical algorithms to alternatively 
solve the Fokker-Planck equation and the macro particle distribution with dynamical weight. 
Finally, Section 6 is dedicated to our summary and conclusions. 
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2 Description of the beam-Beam Interaction 



In order to describe the beam dynamics in an electron positron storage ring, we introduce 
the equations of motion in the following manner. The beam propagation in a reference frame 
attached to the particle orbit is usually described in terms of the canonical conjugate pairs 

fi« = „(*)_ flW^W ; |3« = ^-ff)^£, (2.1) 

9M = S FM + £(„M^_ fl «#) ; (2.2) 

u=x,z \ as Po V PkO ^kO 

where u = (x,z), s is the path length along the particle orbit, and the index k refers to 
either beam (k = 1,2). In equations Q2.1|) and (|2.2j ) the quantity u^ k ' is the actual particle 
displacement from the reference orbit in the plane transversal to the orbit, p^ is the actual 
particle momentum, and E^ is the particle energy. Furthermore, p^ and E k o are the total 
momentum and energy of the synchronous particle, respectively, and is the well-known 
dispersion function. The quantity 

cf<*> = s - 4 k) Rt (2.3) 

is the longitudinal coordinate of a particle from the k-th beam with respect to the syn- 
chronous particle, where is the angular frequency of the synchronous particle and R is 
the mean machine radius. 

It is known that the dynamics of an individual particle is governed by the Langevin 
equations of motion: 

ds dpi k] v ' ds + u 77 ds ' { ' 

ds dr}W u ^ z u ' ds d a(k) ^ r n ' 

where 

F^ = - P ^A k ^ + ^ d ^\ (2.6) 



i + (3-^ + «£ ) W fe) + E^ fc) « (fc) 



(2.7) 
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A k = d\B k \' + JC 2 \B k \ 3/ \ k (s), (2.i 



2r e e 2 55 r e he 3 e 2 

Here a^j is the momentum compaction factor, K^\s) is the local curvature of the reference 
orbit, and = [b^\ Bf~\ Bjfi^J is the magnetic field. The variable £fc(s) is a Gaussian 
random variable with formal properties: 



(Us)) = ; (&00&(s')> = 8(* ~ A (2-10) 
The hamiltonian part in equations (|2.4j) and ( j2.5|) consists of three terms: 



# (fc) = H (k) + ^(fc) + £go ? (2.11) 



where 



F ° = ~ + w|;^7 cos ( v ^ + $fe0 J' (2 ' 12) 



= ^(# 2 + if )2 ) + ^K»^ )2 + (2.13) 



= A fe 5 p ( S )V fe (x( fc \^,ar«; S ). (2.14) 

The parameter K,^ is the so called slip phase coefficient, h k is the harmonic number of the RF 
field and AE k0 is the energy gain per turn. The coefficients G'^s) represent the focusing 
strength of the linear machine lattice, S p (s) is the periodic delta-function, while and 
V k (x^ k \ z^ k \ cr^; are the beam-beam coupling coefficient and the beam-beam potential, 
respectively. The latter are given by the expressions: 

_ r e N 3 ^ k 1 + f3ko/3(3-k)o fn 
M — ^ ' { Z - Lb ) 



V k , z^ k) , a (fc) ; s ) = J dxdzdaQ k (u (k) - u, <r (fc) - a; sj p 3 _ fe (u,a;s), 
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(2.16) 



where N k is the number of particles in the k-th beam and the Green's function Qk(u,a; s) 
for the Poisson equation in the fully 3D case, in the ultra-relativistic 2D case and in the ID 
case can be written respectively as: 



x (k) _ x } _|_ f z {k) _ ^ 



+ (c?(*) - a + 2s 



-1/2 



u (k) - u, cr (fc) - a; s I <; 



5{a {k) -a + 2s) In 



[2.17) 



2n5(a^ -a + 2 s y(z^ - z) 



In what follows we focus on the two-dimensional case, entirely neglecting the longitudinal 
dynamics. Let us write down the Langevin equations of motion ( |2.4j) and ( p.5|) once again 
in the following form: 



ds 



Xis) 



~~ds~ ~ +Fb +Fr ' 



(2.19) 



where 



x (fc) = Mk) ^ 



P w = 



(2.20) 



£! 2 



i? 2 



(2.21) 



is the (external) force acting on particles from the k-th beam, that is due to the linear 
focusing properties of the corresponding confining lattice. Furthermore, 



= X k 5 p (s) 



dV k dV k 



Q2(k) ' Q^(k) 



(2.22) 



is the beam-beam force and 



R 



a (hU) dD i k) Hk) dD^\ 



(2.23) 



is the synchrotron radiation friction force with a stochastic component due to the quantum 
fluctuations of synchrotron radiation [cf expression (|2.8| )1. 
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3 The Nonlinear Sf Method 

It can be checked in a straightforward manner that the Klimontovich microscopic phase 
space density 



i N k 

A(x, p;s) = -^4 x - xi fc )( s )]4 P - pW( a )] (3.1) 



satisfies the following evolution equation: 



^ + P • V,A + (F? + Fg } ) • V p f k + V p ■ (Fg>/ fc ) = 0, (3.2 



9s 



where jx^(s) , p^(s) j is the trajectory of the n-th particle from the k-th beam. Next we 
split the MPSD f k into two parts according to the relation: 



/ fc (x, p; s) = /«)(x, p; s) + <J/ fc (x, p; s), (3.3) 
where f k o is a solution to the equation 



df, 



ds 



- + P • Vs/ko + (F? J + Fg) • V p / fc0 + V p • (FgVfco) = 0. (3.4) 



The quantity F^g in Eq. ( |3.4| ) is the linear part of the beam-beam force F^ . The beam- 
beam force should be calculated with the on-coming beam distribution /(3_mo- In what 
follows it will prove convenient to cast the beam-beam force into the form: 

F<*>=F«+F«+^« (3-5) 

where F^q is the nonlinear (in the transverse coordinates) contribution calculated with 
/(3-fc)0j while 5¥g denotes the part of the beam-beam force due to 5 fa-k- 
it is worthwhile to note here that the representation ( |3.3| ) is unique, embedding the basic 
idea of the 5f method. However, one is completely free to fix the /o part, which usually 
describes those features of the evolution of the system one can solve easily (and preferably in 
explicit form). In the next Section we show that f k0 , averaged over the statistical realizations 
of the process £fc(s) satisfies a Fokker-Planck equation and find its solution. 

Subtract now the two equations (|3.2|) and fl3.4Q to obtain an equation for the 6f k 



96 fk + p • VJf k + (F« + F«) • V p 5f k + V p • (F^5f h 



ds 
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(3.6) 



The next step consists in denning the weight function that is relative to the total distribution 

as 



Substituting 



W fc (x,p; s) 



£/fc(x,p; s ) 
/fc(x,p; s) 



(3.7) 



8fk = W k f k 



fk 



fl 



kO 



1- w fc 



(3i 



into ( |3.6|) and taking into account ( p.2|) we finally arrive at the evolution equation for the 
weights: 



^ + p ■ V x W k + (P« + F« + P«) ■ V p W fc 



fi 



kO 



-(5F 



(k) 

B 



NO 



(3.9) 



Equation (|3.9|) can be solved formally by the method of characteristics. The first couple 
of equations for the characteristics are precisely the equations of motion ( [2.18 ) and ( |2.19|) . 
Suppose their solution (particle's trajectory in phase space) {x(s) , p(s)} is known, and let 
us write down the last one of the equations for the characteristics 



dW h 



W k -l ds 



kO 



— {5F^ + F^o ) ■ Vp/fco 



(*) 



(3.10) 



x,p > trajectory 



Note that its right-hand-side is a function of s only, provided x and p are replaced by par- 
ticle's trajectory in phase space {x(s) , p(s)}. Therefore equation ( pM(J| ) can be integrated 
readily to give: 



W k (s) = 1 + [W k (s ) 



F$(*)] • V p f kQ (a) 



(3.11) 



x(cr) , p(<x) 



8 



4 The Fokker-Planck Equation 



To derive the desired equation let us define the distribution function JF fc0 (x, p; s) and the 
fluctuation 5fko(x., p; s) according to the relations: 



-^fco(x,p; s) = (/fc (x,p; s)) ; <f/fco(x, p; s) = / fe0 (x, p; a) - ^(x, p; s), (4.1) 

where (• • ■) implies statistical average. Neglecting second order terms and correlators in dfko 
and <5/(3_fc)o that generally give rise to collision integrals, we write down the equations for 
F k Q and 6f k0 

<m ° + p ■ V,;F fc o + (F? } + F$) ■ V p ^ + V p ■ (F^JFfco 



5s 



= -V p -(Fg ) e fc (s)5/ fc0 ), (4.2) 

ds 



-V p • (FS^fcCa)^) + 0(6 f M ), (4.3) 



— ffc) ~ (k) 

where F R and F^ denote the deterministic and the stochastic parts of the radiation fric- 
tion force F^ respectively. Moreover, the force F^J should be calculated now with the 
distribution function jF fc0 . Equation (|4.3|) has a trivial solution 



*/*o(s) = -V p • | daFg^a - <r)&( s - a).F fc0 ( s - a), (4.4) 
o 

which is substituted into equation (|4.2| ) yielding the Fokker-Planck equation: 
()s + P • V^o + (F? } + Fg}) ■ Vp^jw + V p ■ (Fg^fco 



dJ^ko 



V p 



F^Vp^F^^o)]. (4.5) 



In order to carry out the 5/ method effectively, it is important to find an equilibrium 
solution of f (or very slowly varying solution) so that the evolution of Sf is separate in time 
scale from that of / . In the following we discuss the equation and the solution of the f 
distribution. 

For the sake of simplicity, in what follows bellow in this Section, we consider one dimen- 
sion only (say x), since the results can be easily generalized to the multidimensional case, 



9 



provided the x-z coupling is neglected. Let us write down the Fokker-Planck equation (|4.5|) 
in the simplified form: 



OS 



F k (s)x— — = T k — (pFko) + V k , 

ox op op op 2 



where 



(4.6) 



2nR 



2ttR 



2ttR 



ds\B k (s) 



ds\B k (s)\ (p kx (s 



(4.7) 



^M + A^( S )AW(, 



/•i(.s) = + A^( S )A^( S ) ; ^(a)* = 



linear part 

Let us seek for a solution of the Fokker-Planck equation ( |4.6| ) in the general form: 



(4.8) 



T k o{x,p\ s) = a fc (s)exp 



j k {s)x 2 + 2a k (s)xp + (3 k (s)p 2 



2e 



xO 



(4.9) 



where e^, is a scaling factor with dimensionality and meaning of emittance. Direct substi- 
tution of (|4.9p into fl4.6|) and equating similar powers (up to second order) in x and p yield 
the following equations for the unknown coefficients: 



da k 
ds 



Vkak ( l ~w)' 



(4.10) 



da k ~ / 2f3 k 

— = F fc /^ fc - 7^ + r fe a fc I 1 



(4.11) 



f =-2* +2 ta(i-^; 



(4.12) 



(is 



'A 



(eg)' 



(4.13) 



where 
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P, 



(eg) 



£>fe 



(4.14) 



is the equilibrium /^-function. 

It is important to note that when the damping vanishes (r& = 0) the above equations are 
exactly the same as the well-known differential equations for the Courant- Snyder parameters. 
In this sense the functions a k , p k and 7^ can be regarded as a generalization of the Courant- 
Snyder parameters in the case when radiation damping and quantum excitation are present. 
The well-known quantity 



Zu = det 



Ik Oik 



(3klk - oi 2 k 



\ oc k (3k / 

is no longer invariant. It is easy to check that its dynamics is governed by the equation 



(4.15) 



ds 



2T k I k [ 1 



Pk 



Pl eq) - 



Comparison between equations (|4.10|) and fl4.16|) shows that 



(4.16) 



ojfc(s) = C k o\ll k (s) 



(4.17) 



with Cfco an arbitrary constant as it should be. Therefore the solution (|4.9|) takes its final 
form 



T k0 (x,p; s) 



Tu(s] 



2lT€ 



X'O 



exp 



Jk(s)x 2 + 2a k (s)xp + P k (s)p 2 

ze x0 



(4.18) 



Let us define now the dimensionless envelope function a k according to the relations 



o~ k 



he 



a k 



P 



Pk 



ke 



P, 



(eg)' 



(4.19) 



Manipulating equations ( 4.11]) , (|4.12|) and ( |4.13| ) for the generalized Courant-Snyder param- 
eters one can eliminate a k and 7^ and obtain a single equation for the envelope a k , which 
combined with equation ( 4.1U| ) comprises a complete set: 



d 2 a k da k 
+ J- fe- 



ds 2 



ds 



F k a k 



o(eq)2 3 
Pk a k°k 



(4.20) 
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dak r / 
— = 1 k ak[l - 



2 2 
a k a k 



(4.21) 



By solving equations ( |4.20|) and (|4.21 ) one can obtain a complete information about the 
evolution of the Tko part of the distribution function. However, solving the above system 
of equations for the beam envelopes and amplitudes of the distributions is not an easy task. 
For that purpose we develop in the next Section a numerical scheme which is more suited 
for direct code implementation. 



5 Numerical Algorithm 

In the previous Sections, we have established the theoretical foundation of the nonlinear 5f 
method for the beam-beam interaction. In this Section we will apply those results to outline 
numerical algorithms suitable for computer simulation. 

Starting with Eq. ( |3.4| ), because the forces in the equation both from lattice and the 
on-coming beam are linear, its solution is well known Gaussian distribution (for example as 
shown in the previous Section in the one- dimensional case) 



1 



T exp 

2 



-z 

2 



T 



(s.i; 



2vrdet (£ fe ) 

where is the matrix of the second moments for the distribution and z is a vector in the six 



dimensional phase space. Based on the method of the beam-envelope |I4" |, the propagation 
of J-"ko can be represented as the iteration of the matrix, 



Mk-tf -Ml + Dk, 



(5.2) 

where A4k is the one-turn matrix including the linear beam-beam force of the on-coming 
beam, and the radiation damping and D k is the one-turn quantum diffusion matrix. Both 
Aik and Dk can be extracted from the lattice using for example the LEGO code [16| . 
However, there is a difference compared to the situation of a single storage ring, namely, 
we have to simultaneously iterate the Gaussian distribution for both beams, since the linear 
map hAk depends on the beam size of the other beam. 

Combining Eqs. (|3TT| ) and (|3.7| ), the perturbative part of the beam distribution 5fk has 



a representation in terms of macro-particles 



=1 



1 



s)S 



P-Pi fe) ( 



where W^ n \s) is the dynamical weight of the n-th particle from the fc-th beam. 



(5.3) 



As a part of the solution for Eq. ( |3.9|) , the propagation of the particle coordinates in 
phase space is the same as the conventional PIC code || provided that the beam-beam force 
is the sum of the two parts from both Tko and Sfk- 
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For the T\& part, we can apply the well known Erskine-Bassetti formula [HJ for a Gaus- 
sian beam. The force due to the Sfk is obtained by solving the two-dimensional Poisson 
equation. In addition to the change of the coordinate, the weight of the particle should be 
propagated according to Eq. ( |3.11|) . The weight should be updated after the change of the 
coordinate since the change of the weight depends on the trajectory of the particle. 



6 Summary 

We have developed an efficacious algorithm for simulating the beam-beam interaction in a 
synchrotron collider with (or without) synchrotron radiation. The nonlinear 5f method has 
been introduced into the evolutionary description of subtle changes of the counter stream- 
ing distribution of the colliding beams over many revolutions. The overall equation that 
describes this evolution is the Fokker-Planck equation (with the radiative process and quan- 
tum fluctuations). In order to isolate the 5f distribution from the average distribution, we 
analyze the solution of the Fokker-Planck equation. Obtained is a form of solution in which 
the time dependence is parameterized through a slow evolution (slow compared with the 
changes in the 5f distribution due to the individual beam-beam interaction) in the Courant- 
Snyder parameters and the emittance of the beam. This algorithm will enhance the analysis 
capability to scrutinize greater details and subtle effects in the beam-beam interaction than 
the PIC version which has been widely deployed ||. 

The current algorithm as well as the previous one || have been developed with an 
immediate application to the PEP-II B-factory collider. The code || has already been 
applied to describe the beam-beam interaction in the PEP-II with unprecedented accuracy 
and reproduction faithfulness, and will be sufficient to study the overall dynamics such as 
the analysis of resonance instabilities and associated luminosity functions. It is anticipated, 
however, that the numerical noise associated with the PIC will require either an inordinate 
amount of macro-particle deployment or a level of noise high enough to mask some minute 
phase space structure that may manifest in subtle but important long-time evolution of the 
beam such as particle diffusion. It is here that the current algorithm will cope with the 
problem. 
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